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Vicente Garzo* and Francisco Vega Reyes^^ 
Departamento de Fisica, Universidad de Extremadura, E-06071 Badajoz, Spain 
(Dated: February 10, 2012) 

A recent segregation criterion [V. Garzo, Phys. Rev. E 78, 020301(R) (2008)] based on the 
thermal diffusion factor A of an intruder in a heated granular gas described by the inelastic Enskog 
equation is revisited. The sign of A provides a criterion for the transition between the Brazil-nut 
effect (BNE) and the reverse Brazil-nut effect (RBNE). The present theory incorporates two extra 
ingredients not accounted for by the previous theoretical attempt. First, the theory is based upon 
the second Sonine approximation to the transport coefficients of the mass flux of intruder. Second, 
the dependence of the temperature ratio (intruder temperature over that of the host granular gas) 
on the solid volume fraction is taken into account in the first and second Sonine approximations. 
In order to check the accuracy of the Sonine approximation considered, the Enskog equation is 
also numerically solved by means of the direct simulation Monte Carlo (DSMC) method to get 
the kinetic diffusion coefficient Do- The comparison between theory and simulation shows that the 
second Sonine approximation to Do yields an improvement over the first Sonine approximation when 
the intruder is lighter than the gas particles in the range of large inelasticity. With respect to the 
form of the phase diagrams for the BNE/RBNE transition, the kinetic theory results for the factor 
A indicate that while the form of these diagrams depends sensitively on the order of the Sonine 
approximation considered when gravity is absent, no significant differences between both Sonine 
solutions appear in the opposite limit (gravity dominates the thermal gradient). In the former case 
(no gravity) , the first Sonine approximation overestimates both the RBNE region and the influence 
of dissipation on thermal diffusion segregation. 
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I. INTRODUCTION 

The understanding of the physical mechanisms in- 
volved in the segregation of an intruder in a granular fluid 
is perhaps one of the most important open challenges of 
granular flow research. This problem has spawned a num- 
ber of important experimental, computational, and the- 
oretical works in the field of granular media [1]. Among 
the different mechanisms proposed to describe (size) seg- 
regation, thermal diffusion becomes the most relevant if 
the system fulfills the conditions of a granular gas. In this 
case, kinetic theory, properly modified to account for the 
inelasticity of collisions, has proven to be a reliable tool 
to analyze the dynamics of the intruder. 

Thermal diffusion (or thermophoresis in its single- 
particle manifestation [2] ) is the transport of matter due 
to the presence of a thermal gradient. As a result of the 
motion of the components of the mixture, a steady state 
can be reached in which the separating effect arising from 
thermal diffusion is balanced by the remixing effect of or- 
dinary diffusion. As a consequence, partial segregation is 
observed and described by the so-called thermal diffusion 
factor A. While this phenomenon has been widely stud- 
ied in ordinary gases and liquids [3], much less is known 
in the case of granular mixtures. It must be noted that in 
the latter case thermal diffusion can appear in vibrated 
systems even in the absence of an external imposed tem- 
perature gradient, as a consequence of inelasticity. In this 



case (energy supplied by vertical walls) , the mean kinetic 
energy of the grains decays away from the vibrating plate 
giving rise to a (granular) temperature gradient. 

In a non-convecting steady state with gradients only 
along the vertical direction (z axis), A is defined by the 
relation [4] 
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where T is the (granular) temperature, and uq and n are 
the number densities of the intruder (or tracer particles 
in a binary mixture) and gas particles, respectively. If 
one assumes that gravity g and the temperature gradi- 
ent point in parallel directions (i.e., the bottom plate is 
hotter than the top plate, InT < 0), then the intruder 
rises with respect to the gas particles if A > while the 
opposite happens if A < 0. When the intruder is larger 
than the gas particles, the former situation is referred to 
as the Brazil-nut effect (BNE) while the latter is called 
the reverse Brazil-nut effect (RBNE). Therefore, the sign 
of the thermal diffusion factor provides a criterion for 
the transition between the BNE and the RBNE by vary- 
ing the parameters of the system (intruder plus granular 
gas). 

A segregation criterion based on the knowledge of the 
factor A has been recently obtained [5] from a solution 
of the (inelastic) Enskog kinetic equation that applies 
to first order in the spatial gradients (Navier-Stokes or- 
der). In contrast to previous theoretical attempts [6, 7], 
the approach is not limited to near elastic particles and 
takes into account the combined effect of thermal gradi- 
ent and gravity. On the other hand, the results reported 
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in Ref. [5] are based on two simplifying assumptions. 
First, although not explicitly stated, they were obtained 
by neglecting the dependence of the temperature ratio 
7 = Tq/T between the intruder and the gas particles on 
the volume fraction 0. This assumption might be more 
questionable as the gas becomes denser (see for instance. 
Fig. 1 below). Second, the results of Rcf. [5] were derived 
by using the first Soninc approximation to estimate the 
transport coefficients associated with the mass flux of in- 
truder. However, recent results [8] obtained in the tracer 
limit for an undriven granular gas clearly show that the 
accuracy of the first Sonine solution can be worsen for 
small values of the coefficients of restitution and/or dis- 
parate values of the mass and size ratios [9] . The ques- 
tion arises then as to whether, and if so to what extent, 
the conclusions drawn from Rcf. [5] may be altered when 
the above two new ingredients (density dependence of 
the temperature ratio and second Sonine correction) are 
accounted for in the theory. In this paper, we address 
this question by determining the thermal diffusion factor 
from the first and second Sonine approximations. 

It is important to remark that the density dependence 
of the temperature ratio affects to both Sonine approxi- 
mations. In this sense, what is new in the present paper 
is not only that the calculations are carried out to higher 
orders, but also that the contributions coming from the 
term dtpj are completely accounted for in the first and 
second Sonine solutions. This fundamental dependence 
was not considered in the previous works on thermal dif- 
fusion for driven [5, 10] and undriven [8] granular gases. 
As the results show, both improvements (density depen- 
dence of temperature ratio and second Sonine correction) 
have a major impact on the physics of the system for 
large mechanical differences between host and intruder 
particles. 

On the other hand, it must be pointed out that the seg- 
regation criterion derived here considers the temperature 
gradient as an input and not created by the inelasticity 
of grains. Thus, if the temperature gradient has a given 
form, then the segregation criterion has a resultant form. 
However, it is known [11-14] for vertically vibrated gran- 
ular systems that after the decrease in the value of the 
granular temperature as a function of height above the 
floor, the temperature profile possesses a minimum above 
which the temperature increases as a function of height. 
Thus, given that dzT < in Eq. (1), our segregation 
criterion can be useful to analyze situations where the 
system is sufficiently small (shallow layers) so that the 
minimum in the temperature profile is not reached or is 
very close to the top of the sample. 

The plan of the paper is as follows. First, the thermal 
diffusion factor A is evaluated in Sec. II by using a hydro- 
dynamic description. This factor is expressed in terms of 
the pressure p of the gas, the transport coefhcients Dq, 
D, and associated with the mass flux of intruder, and 
the (reduced) gravity g* . As for ordinary gases [15], the 
above transport coefficients obey a set of coupled linear 
integral equations that can be approximately solved by 



using a Sonine polynomial expansion. These coefficients 
are determined in Sec. HI by retaining terms up to the 
second Sonine approximation. Some technical details of 
the calculations are relegated to the Appendix A. To 
assess the reliability of the first and second Sonine ap- 
proximations, we compare in Sec. IV the kinetic theory 
predictions for the kinetic diffusion coefficient Dq with 
numerical simulations of the Enskog equation by using 
the DSMC method [16, 17]. As in the undriven case [8], 
the coefficient Dq is computed from the mean-square dis- 
placement of intruders immersed in a heated dense gran- 
ular gas. The knowledge of Dq, D, and D^ allows one 
to express the factor A in terms of the parameter space 
of the system: the mass (rrio/m) and diameter ((To/c) 
ratios, the solid volume fraction (f> and the coefficients of 
restitution a and ao characterizing dissipation of gas-gas 
and intruder-gas collisions, respectively. In order to as- 
sess the impact of the above parameters on segregation, 
in Sec. V the form of the BNE/RBNE phase diagrams 
in the {ctq / a, mo / m}-plane is investigated by varying the 
parameters of the system in the case of hard spheres with 
a common coefficient of restitution (a = ap). Two differ- 
ent limit situations arc mainly analyzed: (i) absence of 
gravity and (ii) thermalized systems (gravity dominates 
the temperature gradient). The paper is closed in Sec. 
VI with a discussion of the results. 



II. THERMAL DIFFUSION OF AN INTRUDER. 
HYDRODYNAMIC DESCRIPTION 

Let us consider a binary mixture of inelastic hard disks 
{d — 2) or spheres {d = 3) where the concentration of 
one of the species (of mass tuq and diameter ctq) is very 
small compared to that of the other (excess component 
of mass m and diameter a < ao). The inelasticity of 
collisions among gas-gas and intruder-gas is accounted 
for by (constant) coefhcients of normal restitution a and 
Qfo, respectively. The system (gas plus intruder) is in 
presence of the gravitational field g = —ge^ , where g is a 
positive constant and is the unit vector in the positive 
direction of the z axis. 

As mentioned in the Introduction, we consider an in- 
homogeneous nonconvecting steady state with only gra- 
dients in the z direction. Since no shearing flows are 
present, then the pressure tensor Pij of the gas is di- 
agonal, namely, Pij ~ pSij, where p is the hydrostatic 
pressure. In this case, the momentum balance equation 
for the gas becomes [10] 



where p = mn is the mass density of the gas particles. 
In the context of the Enskog equation, the hydrostatic 
pressure p is given by [18], 

p = nr[l + 2^^-2^0(1 + a], (3) 

where xi4') is the contact value of the pair cor- 
relation function for the granular gas and (j) = 
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\j^d/2 i2d-ifiY{d/2)]na'^ is the solid volume fraction. Ac- 
cording to the expression (3), the spatial dependence of 
p is through its dependence on the number density n 
(or equivalently, the volume fraction 0) and the granu- 
lar temperature T. Thus, in dimensionless form, Eq. (2) 
yields 
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where p* = p/nT , (3 = d^{(j)p*), and g* = pg/nd^T < is 
a dimensionless parameter measuring the gravity relative 
to the thermal gradient. In addition, since the mean flow 
velocity of the gas vanishes, the mass balance equation 
for the number density no [10] implies jo,2 = Oi where 
jo,z is the mass flux of intruders. 

To close the determination of the thermal diffusion fac- 
tor A, a constitutive equation for the mass flux jo,^ is 
needed. To first order in the spatial gradients (Navier- 
Stokes approximation), the constitutive equation for the 
mass flux jo,z is [18] 
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where Dq is the kinetic diffusion coefficient, D is the mu- 
tual diffusion coefficient, and is the thermal diffusion 
coefficient. The condition — leads to the relation 
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where we have introduced the reduced transport coef- 

)a = {mliy/pT)Do, and 
na'^^^ \J 2T I m is an ef- 



ficients r>^* = {pi^/noT)D'^, D*Q = {mlv/pT)Do, and 
D* = {mav/noT)D. Here, v 
fective collision frequency. 

The explicit form of A can be finally obtained from 
Eqs. (4) and (6). The result is [5] 



A = 



PD'^* - {p* + g*){Dl+D*) 



(7) 



It is quite apparent that in order to assess the impact 
of the parameters of the system (masses, sizes and co- 
efficients of restitution) and the volume fraction <(> on 
the thermal diffusion A, the explicit forms of the diffu- 
sion coefficients D* , Dq and D^* are needed. This can 
be achieved by solving the Enskog kinetic equation by 
means of the Chapman-Enskog method [15]. 



III. ENSKOG KINETIC THEORY. FIRST AND 
SECOND SONINE APPROXIMATIONS TO THE 
DIFFUSION TRANSPORT COEFFICIENTS 

In the tracer limit {no/n — >■ 0), it is expected that the 
state of the granular gas (the solvent) is not affected by 
the presence of tracer particles and that the mutual inter- 
actions of the latter can be neglected as compared with 



their collisions with the particles of the solvent. Conse- 
quently, at a kinetic theory level, the tracer limit implies 
that the velocity distribution function /(r,v,t) of the 
granular gas obeys the (closed) nonlinear Enskog equa- 
tion while the velocity distribution function /o(r, v, t) of 
the intruder obeys the Enskog-Lorentz equation [8]. 

Moreover, as in Rcf. [5], in order to maintain the granu- 
lar medium in a fluidized state, an external energy source 
is coupled to each particle in the form of a thermal bath. 
Here, we consider the situation of energy supply through 
random kicks [19]: the particles of the system are sub- 
mitted between collisions to an uncorrelated white noise. 
This external force is written as !Fi = rrii^, where the 
corresponding stochastic acceleration | is chosen to be 
the same for the intruder and the gas particles [20]. The 
associated forcing term in the Enskog equation is repre- 
sented by a Fokker-Planck collision operator [21] of the 
form — \{C,T/m)d'^/dv'^, where C is the cooling rate asso- 
ciated with the granular temperature T. The generaliza- 
tion of the force to the inhomogeneous case is essentially 
a matter of choice and here, for simplicity, we have as- 
sumed that the stochastic force has the same form as in 
the homogeneous case except that now ^ and T are in 
general functions of space and time. This simple gen- 
eralization has been widely used for ordinary gases in 
shearing problems [22]. It must be emphasized that this 
kind of forcing, which has been shown to be relevant for 
some two-dimensional experimental configurations with a 
rough piston [23] , has been usually employed in computer 
simulations to analyze different problems in the case of 
monodisperse systems [24]. 

The application of the Chapman-Enskog method leads 
to the constitutive equation (5) for the mass flux where 
the transport coefficients D, Dq and are given by 
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As for elastic collisions, the quantities ^ (v), Bq (v) and 
Co (v) are the solutions of a set of coupled linear inte- 
gral equations [see Eqs. (A13)-(A15) of Ref. [10] for the 
driven case]. The standard method consists of approxi- 
mating the above quantities by Maxwellians (at different 
temperatures) times truncated Sonine polynomial expan- 
sions. For simplicity, usually only the lowest Sonine poly- 
nomial (first Sonine approximation) is retained [5, 10, 25] 
and the results obtained from this simple approach agrees 
in general well with Monte Carlo simulations [26]. Ex- 
ceptions to this agreement are extreme mass and size 
ratios in the range of large inelasticity, although these 
discrepancies could be mitigated in part if one considers 
higher-order terms in the Sonine polynomial expansion. 
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This has been recently shown [8] in the undriven case 
(without thermostat) at the level of the transport coefR- 
cient -Dg where the second Sonine approximation yields a 
dramatic improvement (up to 50 %) over the first Sonine 
approximation for some extreme cases. For this reason, 
we evaluate here the complete set of diffusion coefficients 
Dq, D* , and D^* up to the second Sonine approxima- 
tion. The procedure to determine the transport coeffi- 
cients Dq, D and follows similar mathematical steps 
as those made in the undriven case and so we refer the in- 
terested reader to Rcf. [8] for more specific details. Here, 
only the final results are displayed. 

The first Sonine approximations I^o[1]j -C'-^*[1], and 
D*[l] for the (reduced) coefficients Dq, D'^* and D*, re- 
spectively, are given by 
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where 7 = Tq/T is the temperature ratio, Af = niQ/m 
is the mass ratio, lo = aa/a is the size ratio, xo is the 
intruder-gas pair correlation function, fiQ is the chemical 
potential of the intruder and 
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Since granular fluids lack a thermodynamic description, 
the concept of chemical potential appearing in Eq. (13) 
could be questionable. The presence of /ip in our theory 
is justified because we want to recover the results derived 
from the revised Enskog kinetic equation for elastic mix- 
tures [27] (see Appendix C of Ref. [25] for an explanation 
of the choice of some functional derivatives appearing 
in the inelastic Enskog theory). Given that the explicit 
form of the chemical potential must be known to evalu- 
ate the diffusion transport coefficients, for practical pur- 
poses, the expression considered here for /.to is the same 
as the one obtained for ordinary mixtures (a = = 1)- 
Although this derivation requires the use of thermody- 
namic relations that only apply for elastic systems, we 
expect that this approximation could be reliable for not 
too strong values of dissipation. More comparisons with 
computer simulations are needed to support the above 
expectation. 

If we require the stochastic acceleration to be the same 
for both species, it is straightforward to show that he 
temperature ratio 7 fulfills [20, 28] 
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FIG. 1. (color online) Plot of A = —(\>d~f<i> versus the solid 
volume fraction for a hard-sphere gas in the case mo/m = 
8 and ao/a = 2. Three different cases are considered: (a) 
a = Qo = 0.9, (b) a = ao = 0.8, and (c) a = ao = 0.5. 



where in the Gaussian approximation the cooling rates 
C* and Co are 
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The derivative 9^7 appearing in Eq. (13) for D* can be 
obtained by taking the derivative with respect to (j) under 
the condition (15). This yields the relation 
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It must be also remarked that the steady state con- 
dition (15) has been also obtained by considering local 
boundary conditions to the Enskog equation (see the Ap- 
pendix of Ref. [28] ) . More specifically, the boundary con- 
dition considered is a sawtooth vibration of one wall such 
that every particle encountering the wall has a reflected 
speed increased by twice the velocity of the wall in the 
component normal to the wall. In the limit that the 
wall velocity is large compared to the thermal velocities 
of each species (gas particles and intruder), the condi- 
tion (15) is recovered. This equivalence suggests that 
the results derived from the relation (15) can be taken 
as plausible first approximations for qualitative compar- 
isons with experimental results [20]. 

The expression (13) for the diffusion coefficient D*[\\ 
differs from the one derived previously in Ref. [5] by the 
presence of the term drj,^. This contribution was implic- 
itly neglected in those calculations. Equation (13) cor- 
rects this previous approximation. In order to assess the 
effect of this new contribution. Fig. 1 shows the depen- 
dence of A = —(hdj,^ on the solid volume fraction (h for 
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FIG. 2. Temperature ratio 7 = Tq/T versus the (common) 
coefficient of restitution ao = a for a hard-sphere system (d = 
3) in the case ao/a = 2, = 0.2 and two values of the mass 
ratio mo/m: mo/m = 2 (sohd hne and circles) and mo/m = i 
(dashed line and squares) . The lines are the theoretical results 
and the symbols refer to the numerical results obtained from 
the DSMC method. 



three different values of the (common) coefficient of resti- 
tution a = ttQ. We observe that, at a given value of a, A 
increases with increasing cj). In addition, at a given value 
of (j), it is quite apparent that A also increases with in- 
creasing collisional dissipation. Therefore, the influence 
of the density dependence of the temperature ratio on 
thermal diffusion is expected to be important for moder- 
ate densities and strong dissipation. 

The evaluation of the second Soninc approximations 
Z?q[2], Z?^*[2], and D*[2] is more involved. As mentioned 
before, previous calculations carried out in the undriven 
case for the second Sonine approach [8] allows one to eas- 
ily extend these expressions when the gas is heated by 
means of an stochastic thermostat. Their explicit forms 
are provided in the Appendix A. In general, the forms of 
Z?q[2], D'^*\2], and D*[2] have a complex dependence on 
the parameter space of the problem (the mass and size ra- 
tios, the solid fraction and the coefficients of restitution). 
In the elastic limit (a = ao = 1) of a three-dimensional 
gas, Eqs. (A1)-(A3) for the above coefficients agree with 
those previously obtained for a gas mixture of elastic hard 
spheres [29]. Moreover, in the case of mechanically equiv- 
alent particles (toq = to,(Jo = cr, ao = a), as expected, 
one obtains £)^*[2] = and D*q[2] = -D*[2]. Both limit 
cases confirm the self-consistency of the expressions de- 
rived here. 



IV. COMPARISON WITH MONTE CARLO 
SIMULATIONS 

It is important to note that, unless the Sonine expan- 
sion is convergent, the introduction of the second order 
in the Sonine polynomial expansion does not guarantee 
a priori the improvement of the analytical results. In or- 
der to gauge the theoretical predictions one would have 
to compare the latter with computer simulations. This is 
the main goal of this Section where the expressions (11) 
and (Al) of the first and second Sonine approximations, 



respectively, for the kinetic diffusion coefficient will 
be compared with Monte Carlo simulations. 

The diffusion coefficient Do of intruders has been ex- 
tracted by solving numerically the homogeneous Enskog 
equation by means of the DSMC method. As in the un- 
driven case [8, 30], the coefficient Do can be obtained 
from the mean square displacement of the intruder after 
a time interval t as [31] 



|(lr(t)-r(0)l^) = ^. 



(19) 



where ]r(t) — r(0)] is the distance traveled by the intruder 
from t = until time t. Equation (19) is the Einstein 
form of the diffusion coefficient. This relation can be 
used also in Monte Carlo simulations of granular gases 
to measure the diffusion coefficient. In an unbounded 
system like ours, the DSMC method has two steps that 
are repeated in each time iteration. In the first step (free 
streaming stage), the velocity of every particle (in- 
truder and gas particles) is changed to 4- w^, where w.; 
is due to the stochastic force J-i, and is randomly drawn 
from a Gaussian probability distribution (see Ref. [32] for 
more details). The second step accounts for the collisions 
among particles. Because the tracer limit (no/n — >■ 0), 
during our simulations, collisions between intruder parti- 
cles themselves are not considered, and when a collision 
between the intruder and a particle of the gas takes place, 
the post-coUisional velocity obtained from the scattering 
rule is only assigned to the intruder. According to this 
scheme, the numbers of particles have simply a statistical 
meaning, and hence they can be chosen arbitrarily. 

The extension of the DSMC method to study the 
diffusion of intruders in a dense homogeneous granu- 
lar gas requires the changes J[/, /] — )■ xJ[f^ f] and 
MhJ] XoMfoJ]- Here, J[fJ] and Jo[/o,/] refer 
to the (closed) Boltzmann and Boltzmann-Lorentz colli- 
sion operators, respectively [33]. For the DSMC method 
to work appropriately, the time step needs to be small in 
comparison with the microscopic time scale of the prob- 
lem (which is set by the inverse of the collision frequency 
ly) and we also need a sufficiently high number of simu- 
lated particles [16]. We have used in the simulations of 
this paper a time step 5t = 2.5x 10~'^i'~^ and N = 2x 10^ 
simulated particles for each species. 

Before studying the diffusion coefficient Do, it is note- 
worthy to compare the theoretical predictions for the 
temperature ratio 7 with computer simulations. Fig- 
ure 2 shows 7 versus the (common) coefficient of resti- 
tution ao = a for a three-dimensional gas for w = 2, 
(j) = 0.2 and two different values of the mass ratio: M = 2 
and M ~ 1/8. For hard spheres, a good approxima- 
tion for the pair correlation function x is provided by the 
Carnahan-Starling form [34] 



(20) 



(1-0)3' 

while the intruder-gas pair correlation function is given 
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FIG. 3. Reduced kinetic diffusion coefficient Do{a)/ Do{l) as a function of tlie (common) coefficient of restitution a — ag 
for a system of hard splieres with lu — 2 and (p = 0.2. The left panel is for M = 1/8 while the right panel is for M — 2. 
The solid lines correspond to the second Sonine approximation and the dashed lines refer to the first Sonine approximation. 
The symbols are the results obtained from Monte Carlo simulations. Here, Do(l) is the elastic value of the kinetic diffusion 
coefficient consistently obtained in each approximation. 




FIG. 4. Reduced kinetic diffusion coefficient Do{a)/ Do{l) as a function of the (common) coefficient of restitution a — ceo 
for a system of hard spheres with uj — 2 and <j} = 0.2. The left panel is for M — 1/5 while the right panel is for M — 1/2. 
The solid lines correspond to the second Sonine approximation and the dashed lines refer to the first Sonine approximation. 
The symbols are the results obtained from Monte Carlo simulations. Here, 1^0(1) is the elastic value of the kinetic diffusion 
coefficient consistently obtained in each approximation. 



by [35] 

_ 1 o ^ 12 ^^ 
^° 1-0 l + a;(l-</.)2 (1+^^)2(1-0)3- 

(21) 

It is quite apparent the excellent agreement found be- 
tween the theory [obtained from the condition (15)] and 
simulation, showing again the accuracy of the approxi- 
mations (16) and (17) to estimate the temperature ratio. 
As expected [20], the heavier particles carry generically 
more kinetic energy than the lighter ones. Moreover, the 
deviations from the energy equipartition increase as the 
mass differences between intruders and particles of the 
gas increase. As we will show later, in general the effect 
of the temperature differences on the thermal diffusion 
segregation is quite significant. 

Let us consider now the kinetic diffusion coefficient. 
Figures 3 and 4 show the ratio DQ{a) / Dq{1) as a func- 
tion of the (common) coefficient of restitution a ~ aQ 
for w = 2 and </> = 0.2 in the case of spheres {d — 3). 
We have reduced Do{a) with respect to its clastic value 



Dq(1) consistently obtained in each Sonine approxima- 
tion. The solid lines are the theoretical results derived 
from the second Sonine approximation while the dashed 
lines refer to the first Sonine approximation. We observe 
that in general, while the second Sonine approximation 
agrees very well with simulation data, some disagreement 
appears with the first Sonine approximation for strong 
dissipation when the intruder is lighter than gas parti- 
cles. In this case, the first Sonine approximation under- 
estimates the kinetic diffusion coefficient. On the other 
hand, in the opposite case (when the intruder is heavier 
than gas particles) the first and second Sonine approxi- 
mations are practically indistinguishable in the complete 
range of values of a explored and both approaches pro- 
vide a good agreement with Monte Carlo simulations. All 
these results clearly confirm the accuracy of the second 
Sonine approximation for the coefficient Dq, even for low 
values of the coefficient of restitution. Similar conclusions 
were obtained in the undriven case [8] for the diffusion 
coefficient, although the quantitative differences between 
both Sonine solutions were more significant than the ones 
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observed here. 



V. PHASE DIAGRAMS FOR HARD SPHERES 

The exphcit dependence of the thermal diffusion factor 
A on the parameter space of the problem can be obtained 
when one substitutes the expressions of the transport co- 
efHcients D*, D*, and D'^* [Eqs. (11)-(13) for the first 
Sonine approximation and Eqs. (A1)-(A3) for the sec- 
ond Sonine approximation] and the pressure p* (and its 
corresponding derivative /3) into Eq. (7). It is quite evi- 
dent that the influence of the parameters of the mixture 
(masses, diameters, density and coefficients of restitu- 
tion) and the (reduced) gravity on the sign of A is rather 
complicated, given the large number of parameters in- 
volved. In particular, the condition A = provides the 
criterion for the BNE/RBNE transition. Given that the 
results show that both /? and Dq are positive, then, ac- 
cording to Eq. (7), the transition criterion is 

pD^* ^{p*+g*)iD*+D*). (22) 

Note that the thermal diffusion factor is in general a non- 
uniform function since it depends on z through its de- 
pendence on the volume fraction 4'{z) and the thermal 
gradient dzT (in the case that the temperature profile 
is not linear). However, since the expression (7) for A 
has been obtained up to the Navier-Stokes order (first 
order in the spatial gradients), our segregation criterion 
only strictly applies for regions where the density and 
thermal gradients are quite small. Under these condi- 
tions, one can assume that the thermal diffusion factor 
is practically constant, and so the criterion (22) can be 
considered as a global feature of the system. In the case 
of arbitrary spatial gradients, the numerical solution (be- 
yond the Navier-Stokes domain) of the Enskog equation 
via the DSMC method would give better quantitative 
agreement with molecular dynamics simulations or ex- 
periments than the Navier-Stokes results reported here. 

As expected, if the impurities are mechanically equiv- 
alent to the host gas, the system is monodisperse and 
thermal segregation does not occur (A = 0). This is con- 
sistent with our results since in this limit case D^* [2] = 
and Dq[2] ~ -D*[2], so that the condition (22) applies 
for any value of 4> and a. For a dilute gas (0 = 0), the 
first Sonine approximation to Eq. (22) simply yields 

9*h-M) = 0, (23) 

and so no segregation occurs in the absence of gravity. 
This result is consistent with the criterion obtained from 
the Boltzmann equation [36]. However, when the second 
Sonine correction to the diffusion coefficients is retained, 
the results show that thermal segregation appears (A 7^ 
0) for a dilute gas even in the absence of gravity. This is 
illustrated in Fig. 5 showing a non monotonic dependence 
of A on the diameter ratio. In addition, when g* ^ 0, Eq. 
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FIG. 5. (color online) Plot of the thermal diffusion factor A 
obtained from the second Sonine approximation as a function 
of the diameter ratio ao/cr for a dilute {(j> = 0) hard-sphere 
gas in the absence of gravity {g* = 0) when the intruder and 
gas particles have the same mass density (mo/m = (ao/a)^). 
Two different cases are considered: a — ao = 0.7 (solid line) 



— 0.9 (dashed line). 
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FIG. 6. (color online) Phase diagram for BNE/RBNE at 
(j> = 0.1 in the absence of gravity with a — 0.8. Points above 
the curve correspond to A > (BNE) while points below the 
curve correspond to A < (RBNE). The dotted line is the re- 
sult obtained in Ref. [5] , the dashed line is the result obtained 
here from the first Sonine approximation while the solid line 
is the result derived from the second Sonine approximation. 
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FIG. 7. (color online) The same as in Fig. 6 but at (j> = 0.5. 
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(23) leads to the criterion 7 = M . This condition also 
agrees with the Boltzniann results [13] of the undriven 
case. 

To assess the impact of the Sonine approximation on 
segregation, henceforth we only consider the physical case 
of hard spheres {d = 3) with a common coefhcient of 
restitution {a = ap). This reduces the parameter space of 
the problem to five dimensionless parameters: g* , cro/fi 
mo/m, (j) and a. In the case of hard spheres, The expres- 
sion for the chemical potential of the intruder consistent 
with the approximation (21) is [37] 



T 



HnoXl) - ln(l -<P) + 3w-^ 

1-0 



ln(l - (j)) 
21n(l - (f) 



(1-0)2 
0(1-60 + 302) 



1 



, (24) 



where \q{T) is the (constant) de Broglie's thermal wave- 
length [37]. The forms (20), (21), and (24) for x, Xo, 
and /xq, respectively, are the same as those used in our 
previous calculations [5, 10]. 

Figure 6 shows a phase diagram in the (mo/m,CTo/cr) 
plane for a low-density gas (0 = 0.1) in the absence 
of gravity {g* = 0) for a = 0.8. This situation (ther- 
mal gradient dominates over gravity) can be achieved in 
experiments of granular mixtures subjected to horizon- 
tal vibration where the role of gravity can be ignored 
[38]. The corresponding phase diagram obtained from 
the first Sonine approximation by assuming d^^ = [5] 
is also plotted for comparison. We observe first that the 
latter assumption has no significant effect on the first 
Sonine prediction since both results (dotted and dashed 
lines) practically coincide in the range of diameter ratios 
studied. On the other hand, although the first Sonine 
approximation reproduces qualitatively well the trends 
of the phase diagram, the former overestimates dramati- 
cally the predictions of the second Sonine approximation, 
especially at large mass and size ratios. The results also 
show that the quantitative discrepancies between both 
Sonine solutions decrease as the volume fraction in- 
creases. This is illustrated in Fig. 7 for the same system 
as in Fig. 6. The above results suggest that the Sonine ex- 
pansion for thermal diffusion exhibits a poor convergence 
when gravity is absent and the gas is moderately dense 
(say for instance, < 0.2). In particular, the first Sonine 
approximation turns out to be a very poor approximation 
to A in the case of dilute gases. This conclusion agrees 
with the results obtained many years ago by Kincaid et 
al. [4] for ordinary binary mixtures since they concluded 
that the second Sonine approximation is much better ap- 
proximation than the first one. On the other hand, given 
that the first Sonine approximation to Dq agrees well 
with simulation data (see the right panel of Fig. 3) in 
this parameter region (i.e., for mass ratios larger than 
one), the poor convergence of the Sonine approximation 
in this region is mainly due to the coefficients D and . 



We consider now the opposite limit \g* \ — >■ 00, namely, 
when the temperature of the bed is assumed uniform so 
that the segregation of intruder is essentially driven by 
gravity. This situation (gravity dominates over the tem- 
perature gradient) has been previously studied by several 
authors [6, 7] by using kinetic theory and by means of 
computer simulations [39] and experiments [40]. In this 
limit, the condition (22) reduces to Dg + D* = 0. Figure 
8 shows a phase diagram for = 0.2 and a = 0.9 for 
this limit case. It is apparent that, in contrast to the 
case g* — 0, the RBNE regime appears essentially now 
for both large mass ratio and/or small diameter ratio. 
Regarding the influence of the assumption 9^7 = in 
the first Sonine approximation [5] and of the order of the 
Sonine approximation used. Fig. 7 shows clearly that the 
form of the phase diagram is practically independent of 
the approach used since the three curves collapse in a 
common curve. 

Next, we analyze the influence of inelasticity on ther- 
mal diffusion by considering only the most accurate ap- 
proach (the second Sonine approximation). To illustrate 
this effect, two values of a are considered in Fig. ?? for 
= 0.1 and = 0. We see that, in the absence of grav- 
ity, the main effect of collisional dissipation is to reduce 
the size of the BNE, in contrast to what happens for fl- 
nite concentration [41]. In addition, comparison with the 
results obtained for a = 0.5 assuming energy equipar- 
tition (To = T) shows that the impact of nonequiparti- 
tion of granular energy on segregation is important (when 
g* = Q), especially as the size ratio increases. However, 
the influence of nonequipartition is smaller than the one 
previously found from the first Sonine approximation (see 
for instance. Fig. 2 of Ref. [5]). Moreover, a comparison of 
the second Sonine approximation (not shown here) with 
previous theories [6, 7] in the limit \g* \ — > 00 shows simi- 
lar discrepancies as those previously reported (see Fig. 4 
of Ref. [5] ) . This is consistent with the results presented 
in Fig. 6, where the first and second Sonine solutions are 
practically indistinguishable. 



VI. SUMMARY AND DISCUSSION 

In spite of its practical importance, the understanding 
of particle segregation within polydisperse, rapid granu- 
lar fiows is still an open problem. This is due essentially 
to the complexities that arise during the derivation of 
kinetic theory models. Most of the previous theoretical 
works are based on theories that consider systems consti- 
tuted by nearly elastic particles [6, 7, 42, 43], assume an 
equipartition of granular energy [6, 42, 43] and/or they 
are restricted to dilute gases [13, 14, 36]. This paper has 
addressed the problem of segregation by thermal diffu- 
sion of an intruder in a driven granular gas. The analy- 
sis is based on Navier-Stokes constitutive equations with 
diffusion coefficients derived from the revised Enskog ki- 
netic theory [18]. The theory is not restricted to nearly 
elastic spheres, considers the influence of the nonequipar- 
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FIG. 8. (color online) Phase diagram for BNE/RBNE for 
(fi = 0.2 in the absence of a thermal gradient {\g*\ — > oo) with 
a = 0.9. The dotted and dashed lines refer to the results ob- 
tained in Ref. [5] and here, respectively, from the first Sonine 
approximation while the solid line is the result derived from 
the second Sonine approximation. 
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FIG. 9. (color online) Phase diagram for BNE/RBNE at = 
0.1 with = and two different values of a. The solid lines 
are the results from the second Sonine approximation while 
the dashed line is the result derived from the latter approach 
for a — 0.5 but assuming To = T. 



tition of granular energy on segregation, and applies to 
moderate values of the solid volume fraction. In addi- 
tion, in contrast to the previous attempt carried out by 
one of the authors of the present paper [5] , the analysis 
incorporates the density dependence of the temperature 
ratio (9^7) and considers the second Sonine approxima- 
tion (two polynomials in the Sonine polynomial expan- 
sion) to the transport coefficients. These new results are 
the most significant contribution of the present work. In 
this context, this paper complements and extends previ- 
ous papers on segregation of an intruder in driven [5, 10] 
and undriven [8] dense granular gases. 

In the steady state with gradients only along the ver- 
tical direction, the sign of the thermal diffusion factor A 
[defined by Eq. (1)] provides information on the tendency 
of the intruder to move towards the colder (BNE) or hot- 
ter (RBNE) plate. The factor A has been evaluated by 
following two complementary approaches. First, by using 
the momentum balance equation along with the consti- 



tutive equation (5) for the mass flux of the intruder, A is 
expressed in terms of the pressure p (and its derivative 
with respect to the solid volume fraction (f>) and the trans- 
port coefficients Dq, , and D. Then, the form of the 
diffusion transport coefficients has been determined by 
solving the Enskog equation from the Chapman-Enskog 
method up to the second Sonine approximation. This fi- 
nally gives A as a function of the mass and the diameter 
ratios, the volume fraction, the coefficients of restitution 
and the (reduced) gravity g* (a parameter measuring the 
gravity relative to the thermal gradient). 

In order to check the reliability of the first and sec- 
ond Sonine approximations, a comparison with Monte 
Carlo simulations of the Enskog equation for the kinetic 
diffusion coefficient Dq has been performed. As for un- 
driven gases [8], the comparison with simulation data 
shows the superiority of the second Sonine solution over 
the first one specially when the gas particles are heavier 
than the intruder. However, the discrepancies found here 
between both Sonine approximations are less important 
than those previously reported for undriven systems [8]. 

The condition A = [see Eq. (22)] provides the seg- 
regation criterion for the transition BNE <^ RBNE. A 
systematic study of the form of the BNE/RBNE phase 
diagrams in the mass and diameter ratio plane has been 
carried out in Sec. V for hard spheres {d ~ 3) in the case 
of a (common) coefficient of restitution a = ao- Regard- 
ing the form of the phase diagrams, we observe first that 
the impact of the term 9^7 on the first Sonine approxi- 
mation is very small (see, for instance. Figs. 6 and 7) so 
that the conclusions made in Ref. [5] are not practically 
altered by the presence of this new term. With respect 
to the effect of the Sonine approximation considered, the 
results reported here show that in general the influence of 
the order of the Sonine solution is much more significant 
in the absence of gravity {g* = 0) than in the opposite 
limit (|g*| — > oo), where both approximations yield prac- 
tically the same diagrams (see, for instance. Fig. 6). In 
particular, when g* — 0, the first Sonine approach clearly 
overestimates the RBNE region (see Fig. 6), while it pre- 
dicts a bigger influence of coUisional dissipation on the 
phase diagram than the one obtained from the second 
Sonine correction (compare Fig. 2 of Ref. [5] with Fig. 
7 of the present paper). Although the accuracy of the 
second Sonine approximation to the diffusion coefficients 
has been only tested in the case of the coefficient Dq, 
it can be reasonably expected that the segregation cri- 
terion derived here in the second Sonine approximation 
for a heated gas compares better with computer sim- 
ulations than segregation criteria from previous works 
in the Navier-Stokes domain (small gradients of density 
and temperature). A previous comparison [13] for dilute 
gases confirms this expectation. Given that the present 
results apply to moderate densities, it is hoped that this 
paper could stimulate the performance of such simula- 
tions. 

An important question addressed partially in this pa- 
per is about the convergence of the Sonine polynomial 
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expansions considered here. This is a quite difficult ques- 
tion, specially in the case of granular fluids where the 
studies of the impact of higher-order terms on transport 
are more scarce than for ordinary gases. In this latter 
case, for instance the analysis of transport properties of 
dense binary mixtures with one tracer component [44] 
indicates that the convergence of the Soninc expansion 
improves significantly with increasing values of the mass 
ratio M . Similar trends have been found here at the level 
of the tracer diffusion coefficient Dq. Unfortunately, the 
lack of available simulation data for the remaining dif- 
fusion transport coefficients D and D"^ in the driven or 
undriven cases prevent us to assess the reliability of the 
first and second Sonine approximations. Since those coef- 
ficients are also involved in the expression of the thermal 
diffusion factor A, no definitive conclusions on the accu- 
racy of their second Sonine forms can be drawn. How- 
ever, the results displayed in Section V for thermal diffu- 
sion segregation seem to indicate that while the conver- 
gence of the Sonine expansion for A is quite good when 
I — )■ oo (see for instance Fig. 8 where both Sonine ap- 
proaches lead practically to the same phase diagrams), 
it does not happen the same for a low-density gas in the 
opposite limit {g* = 0) since there is an abrupt change 
in the form of the phase diagram for BNE/RBNE (see 
for instance Fig. 6). In this situation one should perhaps 
consider higher-order polynomial terms (even beyond the 
second Sonine approximation) or one should consider al- 
ternative analytical routes, such as the so-called modified 
Sonine approximation [45]. This latter method is based 
on a modified version of the first Sonine approximation 
which replaces the Gaussian distribution weight function 
(used in the standard Sonine method) by the homoge- 
neous cooling state distribution. In any case, more com- 
parisons between segregation results derived in the driven 
case between computer simulations and the different ap- 
proximate theories are needed before quantitative con- 
clusions can be offered on the reliability of those kinetic 
theories. 



ACKNOWLEDGMENTS 

The present work has been supported by the Ministe- 
rio de Educacion y Ciencia (Spain) through grants No. 
FIS2010-16587 (VG and FV) and No. MAT2009-14351- 
C02-02 (FV), partially financed by FEDER funds and by 
the Junta de Extremadura (Spain) through Grant No. 
GR10158. 



played in this Appendix. They are given by 



D*[2\ 
where 



(Al) 



, (A2) 



(A3) 



XI = - {Mp* - 7) + -(1 + iofMoXoH'^ + ao), (A4) 



1 A//^ 

• "-(l+c^)V'xo<^(l + ao) 



2{d + 2) M 



X 1^ [(rf + 2){Ml - 1) + (2d - 5 - 9ao)MoM 
+ (d - 1 + 3ao + 6al)M^] + GM^il + aof} , 

(A5) 



2(1 + M) T 



T,no 



(A6) 



1 



^ d //io\ 



X { [(d + %)Ml + (7 + 2d - 9ao)MoM + {2 + d 

+ [{d + 2.)Ml + (2d - 5 - 9ao)-MoA^ + (d - 1 
+3ao + 6a^)A^2j 5)2 _ ^ 2)6'(1 + 9)]. (A7) 

Here, JV[ = m/{m + TOq), AAq = mQ/{m + TOq) ^-nd 
9 = nioT/uiTQ is the mean-square velocity of the gas par- 
ticles relative to that of the intruder particle. In addition, 
when the gas is driven by the stochastic thermostat, the 
coefficients a* and c* are given by [46] 

a* = (^1 + 3^^^0x(l + «)'(2a - 1)) , (A8) 



c* = -3zyr'T?^<^(2x + 0a^x)a(l-a'), (A9) 



(d + 2) 



Appendix A: Second Sonine expressions for the 
diffusion coefficients 



The explicit expressions of the second Sonine approx- 
imations to the diffusion transport coefficients are dis- 



where 



, _ 8 7r('^-i)/2 
d(d + 2)V2r(d/2)^^ 
'd- 1 3 



+ — (d + 8)(l-a) . (AlO) 



2 16 
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Finally, the expressions of the (reduced) collision frcqucn- where G and D arc given by Eqs. (C7) and (C8), respec- 
cies , and can be found in the Appendix C tivcly, of Ref. [8] . 

of Ref. [8] while and vl are given by [47] 

_(d-l)/2 / 1 I , ,\ d.-l 



(All) 



1^ 
3/2 



d(d + 2)r(|) V 2 J ^""Mo 



[C+{d + 2){l + 9)D], (A12) 
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